
//High log
/*
	s = X(i,j)*2^m;
	R(i,j) = pi / (2*agm(1,4/s,100)) - m * log(2);
agm:
	a = 0.5*(x+y);
	g = sqrt(x*y);

	for i=1:it
	    at = a;
	    a = 0.5 * (a+g);
	    g = sqrt(at*g);
	end; 
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define ITER 100
#define AGM 20
#define PI 3.14159
#define LN2 0.69314
#define MBITS 8

int main(){

	float f;
	int i,j;
	float a;
	float g;
	float at;
	float y;
	float lnf;
	float s;
	float x;


	for(i=0; i< ITER; i++){
		f = 1/(1-((float)rand()/RAND_MAX));
		
		s = f*pow(2,MBITS);		
		y = 4.0/s;
		x = 1;
		//agm start
		a = 0.5*(x+y);
		g = sqrt(x*y);
		for(j=0; j<AGM; j++){
			at = a;
			a = 0.5*(a+g);
			g = sqrt(g*at);
		}
		//agm end
		
		lnf = PI / (2*a) - MBITS * LN2;

		printf("%f %f %f %f\n", f, lnf, logf(f),fabs(lnf-logf(f)));
		
	}

	return 0;
}
